;-----------------------------------------------------------------
; 
; bivariate PDF distribution of RN222, SO2, and T  
;
;   - Using a user-defined function for cleaner code
;   - Subselecting an array using its coordinate arrays
;   - Using the "where" function
;   - Using "raster" contour mode for better delination of data
;   - Using "plt_pdfxy"
;-----------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
load "~/wrk/lib/lib_radon.ncl" 


begin 

  datapath = "./" 
  plotpath = "./" 
  plotform = "eps"     

  dist = "RN222_SO2" 
  dist = "RN222_T_P" 
  dist = "RN222_T" 

  sitename = "Freiburg" 
  sitename = "Beijing"
  sitename = "Schauinsland" 
  sitename = "Socorro" 
  sitename = "EUR" 
  sitename = "ASIA2" 
  sitename = "USA" 
  sitename = "AUS" 

  plotname = dist + "_" + sitename 

  flnma = "RA008_200001.nc" 
  flnmb = "RA010_200001.nc" 
  flnmc = "RA007_200001.nc" 
  flnmd = "SO2_200001.nc" 
  flnme = "st_200001.nc" 
  flnmo = "PDF_"+sitename+".nc" 


  if(True) then 

  fla = addfile(datapath+flnma, "r")
  flb = addfile(datapath+flnmb, "r")
  flc = addfile(datapath+flnmc, "r")
  fld = addfile(datapath+flnmd, "r")
  fle = addfile(datapath+flnme, "r")

  rna = fla->RN222
  rnb = flb->RN222
  rnc = flc->RN222
  so2 = fld->SO2
  st  = fle->st

  rna = where(rna.gt.0.,rna,0.) 
  rna = RN222_ionprod(RN222_mBqm3(rna))

  rnb = where(rnb.gt.0.,rnb,0.) 
  rnb = RN222_ionprod(RN222_mBqm3(rnb))

  rnc = where(rnc.gt.0.,rnc,0.) 
  rnc = RN222_ionprod(RN222_mBqm3(rnc))

;-----------------------------------------------------------------
; Off the coast of S. America at low altitudes (Liquid):  
;     15-30S, 75-90W, 1000-700 hPa
;-----------------------------------------------------------------

  latS   =  40. ;; EUR 
  latN   =  60.
  lonL   =  10.
  lonR   =  40.

  latS   =  20. ;; ASIA2 
  latN   =  60.
  lonL   =  60.
  lonR   = 120.

  latS   =  20. ;; USA 
  latN   =  60.
  lonL   = 240.
  lonR   = 300.

  latS   =  20. ;; AUS 
  latN   =  60.
  lonL   = 240.
  lonR   = 300.


  latP   =  40. ;; beijing 
  lonP   = 116. 

  latP   =  48. ;; freiburg 
  lonP   = 7.7  

  latP   =  34. ;; socorro
  lonP   = -106.7 + 360. 

  printVarSummary(rna) 
  printVarSummary(rnb) 
  printVarSummary(rnc) 

;-----------------------------------------------------------------
; Bivariate (Joint) PDFs
; Use coordinate subscripting to select data
;-----------------------------------------------------------------

  klev = 30

  if(sitename.eq."Schauinsland") then 
     klev = 28
  end if 

  if(dist.eq."RN222_SO2") then 

  pdfa = pdfxy(rna(1:,klev,{latS:latN},{lonL:lonR}) \
               ,so2(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  pdfb = pdfxy(rnb(:,klev,{latS:latN},{lonL:lonR}) \
              ,so2(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  pdfc = pdfxy(rnc(:,klev,{latS:latN},{lonL:lonR}) \
              ,so2(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  end if 


  if(dist.eq."RN222_T") then

  pdfa = pdfxy(rna(1:,klev,{latS:latN},{lonL:lonR}) \
               , st(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  pdfb = pdfxy(rnb(:,klev,{latS:latN},{lonL:lonR}) \
              , st(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  pdfc = pdfxy(rnc(:,klev,{latS:latN},{lonL:lonR}) \
              , st(:,klev,{latS:latN},{lonL:lonR}) ,0,0,False)

  end if


  if(dist.eq."RN222_T_P") then

  pdfa = pdfxy(rna(1:,klev,{latP},{lonP}) \
               , st(:,klev,{latP},{lonP}) ,0,0,False)

  pdfb = pdfxy(rnb(:,klev,{latP},{lonP}) \
              , st(:,klev,{latP},{lonP}) ,0,0,False)

  pdfc = pdfxy(rnc(:,klev,{latP},{lonP}) \
              , st(:,klev,{latP},{lonP}) ,0,0,False)

  end if

;-----------------------------------------------------------------
; Set all pdf2=0.0 to _FillValue 
; Better illustrates actual data distribution
;-----------------------------------------------------------------

  pdfa = where(pdfa.eq.0, pdfa@_FillValue, pdfa)
  pdfb = where(pdfb.eq.0, pdfb@_FillValue, pdfb)
  pdfc = where(pdfc.eq.0, pdfc@_FillValue, pdfc)

  system("rm "+flnmo) 

  flo = addfile(flnmo,"c") 
  flo->PDFA=pdfa 
  flo->PDFB=pdfb 
  flo->PDFC=pdfc

  else

  flo = addfile(flnmo,"r") 

  pdfa = flo->PDFA
  pdfb = flo->PDFB
  pdfc = flo->PDFC

  end if 


;-----------------------------------------------------------------
; plots
;-----------------------------------------------------------------

  colormap = "testcmap" 
  colormap = "amwg" 

  plot = new(3,"graphic")

  wks = gsn_open_wks(plotform, plotname)

  gsn_define_colormap(wks,colormap)

  res                      = True
  res@gsnDraw              = False    ; do not draw picture
  res@gsnFrame             = False    ; do not advance frame
  res@cnFillOn             = True     ; turn on color fill
  res@cnFillMode           = "RasterFill"       ; Raster Mode
  res@cnLinesOn            = False    ; no contour lines
  res@cnLineLabelsOn       = False    ; no contour line labels
  res@cnInfoLabelOn        = False
  res@gsnSpreadColors      = True     ; use full colormap
 ;res@cnLabelBarEndStyle   = "ExcludeOuterBoxes"

 ;res@lbLabelBarOn         = False
  res@lbLabelAutoStride    = True

  res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiXAxisOffsetYF      = 0.0 


  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  0.2               ; set min contour level
  res@cnMaxLevelValF       =  2.                ; set max contour level
  res@cnLevelSpacingF      =  0.2               ; set contour spacing
  ;res@cnMaxLevelValF       =  1.                ; set max contour level
  ;res@cnLevelSpacingF      =  0.1               ; set contour spacing

  res@gsnLeftString        = ""
  res@gsnCenterString      = "WCRP1995 "
  res@tiYAxisString        = "Temperature (K)" 
  plot(0) = gsn_csm_contour (wks, pdfa, res)

  res@gsnLeftString        = ""
  res@gsnCenterString      = "SW1998 / 1.6 "
  res@tiYAxisString        = ""
  plot(1) = gsn_csm_contour (wks, pdfb, res)

  res@gsnLeftString        = ""
  res@gsnCenterString      = "MERGE " 
  res@tiYAxisString        = ""
  plot(2) = gsn_csm_contour (wks, pdfc, res)

  resP                     = True                ; modify the panel plot
 ;resP@txString            = title
  resP@gsnPanelRowSpec     = True                   ; tell panel what order to plt
  ;resP@gsnPanelLabelBar    = True
  ;resP@lbLabelAutoStride   = True
  ;resP@lbOrientation       = "vertical"          ; vertical label bar
  gsn_panel(wks,plot,(/3,1/),resP)

;-----------------------------------------------------------------
; Use "plt_pdfxy" to plot
;-----------------------------------------------------------------

;  res@gsnCenterString      = "PDF"
;  plot(0) = plt_pdfxy (wks, pdfa, res)
;
;  resP                     = True 
;  resP@txString            = "" 
;  resP@gsnPanelRowSpec     = True              
;  gsn_panel(wks,plot,(/1,1/),resP)

end                 

